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Judging Model Reduction of Chaotic Systems via Optimal Shadowing Criteria 
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A common goal in the study of high dimensional and complex system is to model the system by 
a low order representation. In this letter we propose a general approach for assessing the quality of 
a reduced order model for high dimensional chaotic systems. The key of this approach is the use of 
optimal shadowing, combined with dimensionality reduction techniques. Rather than quantify the 
quality of a model based on the quality of predictions, which can be irrelevant for chaotic systems 
since even excellent models can do poorly, we suggest that a good model should allow shadowing by 
modeled data for long times; this principle leads directly to an optimal shadowing criterion of model 
reduction. This approach overcomes the usual difficulties encountered by traditional methods which 
either compare systems of the same size by normed- distance in the functional space, or measure how 
close an orbit generated by a model is to the observed data. Examples include interval arithmetic 
computations to validate the optimal shadowing. 

PACS numbers: 05.45.-a 05.10.-a 05.45.Xt 89.75.Hc 
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Model reduction is an important concept found across 
science and engineering. Approximating gross scale fea- 
tures of high dimensional systems is a fundamental ques- 
tion which occupies a great deal of time and energy in 
the study of such disparate mathematical fields as PDE 
theory, time-delay systems, networked dynamical sys- 
tems, and where-ever high dimensional problems natu- 
rally arise from the underlying science from which come 
the models. The POD method for example [l| is a pop- 
ular way to produce a basis set for high-dimemensional 
data from solutons of PDEs, onto which the resulting 
Galerkin projections are optimal in the sense of a fastest 
decaying time-average power spectrum. Underlying such 
techniques there is usually the common thread of min- 
imization of the I2 distance in the functional space be- 
tween the actual system and its reduced order model - 
models are considered best in the Banach space. How- 
ever, for chaotic systems, use of the £2 minimization crite- 
ria to compare the two functions for determining whether 
a model is good may not be relevant, since two functions 
can be close in an underlying Banach space, but exhibit 
dramatically different dynamical properties Q. 

Likewise, a reasonable model, even a perfect model, 
may quickly produce quickly and dramatically differ- 
ent simulation results - it is well known that comparing 
time-series from simulations is an unworkable criterion 
of model comparison due to sensitive dependence. When 
random noise or modeling error is introduced, as is ar- 
guably always the case in practice, even a seemingly per- 
fect model would suffer from conflicting judgements. The 
sensitivity to perturbations prevents us from the compar- 
ing chaotic systems by direct comparison of their trajec- 
tories, since even (almost) identical systems would fail 
such a measure of comparison. See Fig. [1] as an example. 

To judge a model reduction, it is too much to hope that 
a model will be capable to reproduce trajectories of the 
full system, due to the chaotic nature of the system, as 



well as technical details of comparing trajectories which 
arise from systems of different dimensionality We assert 
that such comparisons are meaningless good or bad be- 
cause the expectation that the results will always be bad. 
Instead, we will judge a model to be a good representa- 
tion if its trajectories can numerically shadow trajectories 
of the full system. In this sense, the model is producing 
plausible solutions, if not the actual simulations. 

Shadowing was introduced initially to rigorously verify 
the existence of a true orbit from a model to a computer 
generated orbit which is usually noisy Given a 

noisy orbit p = {pt} 7 a model generated orbit x = {x t } 
is said to e— shadow p if ||x — p||oo = su Pt \ \ x t — Pt\U < 
e From now on these subscripts of norms will be 

omitted unless otherwise specified. 

In terms of judging the model quality, we wish to asso- 
ciate the capability of the model to shadow observation 
with its quality. However, most shadowing techniques 
were developed only to find an arbitrary e— shadowing 
orbit, which may be far from optimal (there may be an- 
other shadowing orbit with a much smaller e), preventing 
us from a good judgement of the model. To overcome this 
ambiguity, we ask: what is the best orbit the model can 
produce, to match the observed orbit? This amounts to 
judge the quality of a model / : D — > D for given obser- 
vation p by the optimal shadowing distance: 
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where D C K. m and x is an orbit of x\ under / [1 91 ] . 
For deterministic systems, this question is equivalent to 
finding the initial point which leads to a true orbit that 
can step-by-step match the noisy orbit best. 

Based on the concept of optimal shadowing, we focus 
on the question of how to understand the quality of a 
model reduction, meaning how well does a model of lower 
dimensionality represent the dynamics of the full system. 
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FIG. 1: Illustration of the difficulty in judging a model by 
comparing orbits directly. In the upper panel, a noisy numer- 
ical orbit p = {ptYtL\ of the logistic map is shown (in black 
square). This orbit satisfies pt+i = 4pt(l —pt) + St where 8t is 
uniformly distributed in [— 2 10 ,2 10 ]. Blue triangle represents 
a true, noiseless orbit z = {z t }tZi with z\ = pi = 0.8724.... 
Note that although z is close to p for initial times, after about 
10 steps they start to diverge. On the other hand, starting 
with si = 0.8723..., we found a true orbit s = {st}f£ 1; shown 
in (red) crosses, which is able to match the entire noisy orbit 
p. Although generated by the same model, p and s appar- 
ently leads to different conclusions about the model quality if 
we were to judge a model by comparing time series naively. 



Since a high dimensional system and its reduced order 
model necessarily generate time series of different dimen- 
sions, there is currently no direct way of comparing two 
such models. Our approach to solve this problem can be 
illustrated by the diagram in Fig. [5] Given a high dimen- 
sional system and its candidate reduced order model: we 
first generate time series from the original system; next, 
dimensionality reduction is performed to extract a low 
dimensional representation of the time series; finally, we 
look for an optimal shadowing orbit from the reduced 
order model to match the low dimensional time series. 
The reduced order model, being a simplification of the 
original one, suffers from two types of inexactness. The 
first type comes from dimensionality reduction, which ac- 
counts for the loss of information in simplifying the ob- 
servation; the second type comes from shadowing, and 
is crucial for assessing the model quality of chaotic sys- 
tems, which here accounts for the capability of the given 
model to generate one orbit that matches the observed 
(low dimensional) time series. 

This approach allows us to quantify the quality of a 
model reduction even for chaotic systems, which is not 
likely to be achieved by traditional methods. Further- 
more, the flexibility in emphasizing in between dimen- 
sionality reduction and shadowing errors allows one to 
adjust the measure of model quality in different situa- 
tions depending on specific applications. 

To illustrate this perspective, we consider the prob- 
lem of modeling a system of coupled chaotic oscillators. 
Coupled oscillators have been studied extensively as pro- 
totypical of complex systems [1, Q , with promising ap- 
plications ranging widely from the modeling of flocking 
behavior [13], to mathematical epidemiology where col- 
lective behavior leads to mean field model of disease dy- 
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FIG. 2: General model reduction design cycle. (Down) A 
large scale system gives rise to a many variate time-series. 
(Bottom Across) Averaging across scales gives rise to a lower 
dimensional system with correspondingly fewer measurable 
variables in the output time-series. (Up) The step of judging 
model quality is usually overlooked in the design cycle. To in- 
fer model quality, in some way the model must be required to 
remind of the full system. Here we advocate that prediction 
is inappropriate due to sensitive dependence to initial condi- 
tions. Instead we suggest an optimal shadowing criterion. 



namics to mention a few. In any of these settings 
where many coupled oscillators may arise, it is natural 
to average across spatial scales so that a model with just 
a few oscillators may be meant to represent the system, 
in the sense that an element of the model may represent 
many elements of the whole. In the much the same way 
as a community analysis of complex networks where the 
topology allows partitioning into groups [l2, 13 1, in dy- 
namical systems we assert that groups of oscillators may 
exist with similar behavior. When a system is modeled by 
a large collection of coupled oscillators the natural ques- 
tion is how might the simplified low order model captures 
similar properties of the original high order system? 

We choose to illustrate our approach by a system of 
coupled quadratic maps, described by: 
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where {z^}i=i,...,n represent a set of coupled oscillators, 
G 5R d is the state of oscillator i at time i; each 
individual oscillator is driven by a discrete logistic dy- 
namics with parameter a^ l \ which allows possible mis- 
match of parameters between different individual oscil- 
lators, which is usually the case for a physical setting; 
the second term describes the effective coupling between 
different oscillators through a discrete Laplacian matrix 
L = [Uj] n xn, where for each i, Y^j=i % = 0; and a is the 
coupling strength. The coupling function has been cho- 
sen to have the same form of the individual dynamics, 
which corresponds to the situation where each oscillator 
receives a direct signal from the output of its neighbors. 
For this high (n x d) dimensional coupled system, sev- 
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eral questions are of particular interest, as initial explo- 
ration for the general problem, and will be answered in 
this letter. Fig. [3] serves as an illustration. 

1 In what sense can we model a coupled identical oscil- 
lator network by a single oscillator? 

2 In what sense can we model a coupled non-identical 
oscillator network by a single oscillator? 

3 In what sense can we model a nearly synchronized clus- 
ter by a single oscillator? 



For question 1, a general criteria is whether the system 
synchronizes or not. When the oscillators synchronize, 
lim t _ i . 00 | — x\ || — > for Vi,j. After transient, all the 
oscillators evolve in the same way, and the second term in 
Eq. ([2]) disappears (there will be no error in dimension- 
ality reduction or shadowing). Any single oscillator i is 



governed by the same dynamics: xft-y 
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Thus we can perfectly model the coupled system by a 
single, low dimensional system: st+i = as t (l — Sj). 
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FIG. 3: Illustration of a model reduction of the coupled oscil- 
lator network. In the first case (top ellipse), all the oscillators 
are the same; in the second case (middle ellipse) the oscillators 
are mismatched; while in the third case, the network consists 
of a cluster of identical oscillators with a few outliers. The 
rectangles represent individual oscillators (the width and color 
are used to highlight the difference of individual oscillators). 
In all cases, we are interested in whether model reduction is 
possible from the original system to a low dimensional system 
(represented by a single oscillator on the right). 

Questions 2 and 3 are intriguing. In these cases, the 
oscillators are unable to completely synchronize, thus a 
single oscillator model may not exactly represent the true 
collective behavior of the coupled system. In particular, 
if one chooses the average trajectory as a low dimensional 
representation of the high dimensional time series, then 
this average variable is governed by 

n n 

- at+1 = Lj2* (i) 4 ) (i-4 i) )-- E h ja (H j) (i-4 j) ), 



which depends essentially on every single oscillator, im- 
plying that the dimension of the system is as high 
as the original coupled system. Even in the situa- 



tion where the oscillators are nearly synchronized 14| : 
limsup t ||a4 — ^tll ~ 0) if one were to use mean-field 
approximation, replacing xf with x t and with a, re- 
sulting in a model: 



a t +i = ax t {l - x t ), 



(4) 



(3) 



then at each step this model generates error (compar- 
ing to the actual average state) which comes from the 
heterogeneity of the individual dynamics, and its effect 
might be tremendous depending on how the heterogene- 
ity distributes among the oscillators. Nevertheless, our 
approach overcomes the difficulty and provide a quanti- 
tative measure of the reduced order model. 

We shall illustrate this for case 2 by the use of op- 
timal shadowing for the average trajectory. As a mat- 
ter of example, we will construct a network of logistic 
oscillators whose individual parameters are drawn 
uniformly from [3.9998,4] in order to emphasize oscilla- 
tor mismatch. We couple those oscillators through an 
Erdos-Renyi network [15[ of n = 1000 and p = 0.1 (the 
probability that any two nodes are joined by an edge), 
with coupling strength a = 0.0075. 

The dependence of optimal shadowing distance de- 
pends upon the parameter a for a one-parameter family 
of reduced models f{x) = ax(l — x). Here we use a finite 
trajectory of length T = 1000 after transients. e opt are 
calculated by use of interval arithmetic, with the excel- 
lent package "INTLAB" [lf|, in order to validate that 
we are representing reasonable upper bounds of the ac- 
tual optimal shadowing distances. Results are shown in 
Fig. [31 for a typical trajectory generated by the original 
network. It is interesting to note the difference between 
using the shadowing criteria in contrast to the usual £2 
criteria: while the model error seems to depend symmet- 
rically on a under the £2 criteria, shadowing is able to 
capture the asymmetry which seems to be more reason- 
able because of the increase of topological entropy for 
increasing a. Shadowing also has the advantage to judge 
how long the reduced order model is valid for the orig- 
inal system (the optimal shadowing distance increases 
non-smoothly when we take longer trajectories), another 
perspective the £2 criteria does not provide. We have also 
obtained similar results in the case of modeling a nearly 
synchronized cluster (case 3), which will be reported in 
a more detailed paper. 

The above example demonstrates the judging of a 
model reduction by measurement of the optimal shad- 
owing distance from a model to the average trajectory 
from the original system. Our choice to use such an av- 
erage was selected to minimize the square distance to 
all other individual trajectories, i.e., the dimensional- 
ity reduction error. To illustrate this perspective, con- 
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FIG. 4: (Color online) Optimal shadowing distances of re- 
duced order model for a coupled oscillator network. This 
system consists of mismatched logistic oscillators coupled 
through a network of n = 1000 nodes and p = 0.1, with cou- 
pling strength a — 0.0075. Blue square corresponds to the 
optimal shadowing distance for a one-parameter family of re- 
duced order models f(x) — ax(l — x); black star corresponds 
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to the £ 2 norm \JJ2 t =i \f( x t) ~ x t +i\ 2 /T. 



FIG. 5: (Color online) Interplay between dimensionality re- 
duction error and shadowing error. In each panel the hori- 
zontal axis corresponds to A and vertical axis corresponds to 
the model reduction error: ("(A) = (1 — n)r](X) + e op t(A), for 
fixed /j,. 



sider a toy example where we have two logistic oscil- 
lators with parameters 3.9998 and 4. coupled through 
a network of Laplacian matrix L — [2,— 2;— 1,1] with 
coupling strength a = 0.25. The dimensionality reduc- 
tion of the time series [x^^x' 2 )] can be represented by 
a convex sum: x = (1 — A)xW + Ax' 2 ). For given A, 
the dimensionality reduction error can be defined as: 

77(A) = H x - x(i) ll 2 /( 2T ) where T is the len S th 

of x. In Fig. [5] we show how one would obtain dif- 
ferent dependence of the model reduction error C(A) 
: 1 — /i)r/(A) + e op t(A) on A. It is interesting to note espe- 
cially in the last panel (lower right of Fig. [5]) that when 
we emphasize purely on the modclability of the low di- 
mensional system, then the trajectory from the single 
oscillator x' 2 ' would induce the best model (among the 
family of models f(x) = ax{l — x)). On the other hand, 
for other choice of /x, the optimal A would change, not 
necessarily equals 1/2, as expected. 

In general it will be interesting to ask such questions as 
in a large network, how shall we take the weighted aver- 
age of individual trajectories to reach an optimal balance 
between dimensionality reduction and shadowing; or how 
nonlinear dimensionality reduction can be adopted in the 
case of generalized synchronization 17J . Some of the re- 
sults will be reported in a future paper. 

To summarize, we have proposed a general approach 
for assessing the quality of reduced order models for high 
dimensional chaotic systems. The key in this approach is 
the unusual application of concepts from shadowing, to- 
ward the optimal shadowing criterion, combined with di- 
mensionality reduction techniques. This approach over- 
comes perhaps overlooked problems inherent with tradi- 
tional methods of comparison which may cither attempt 



to compare systems of the same size by measuring the 
distance in the functional space, or alternatively to mea- 
sure how close an orbit generated by a model is to the 
observed data. Both of these perspectives have funda- 
mental flaws which our optimal shadowing based cost 
function overcomes. 
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